Stochastic Block Models

Network Analysis

Graph Theory

Community Detection

A generative model for random graphs used to find communities of nodes with similar connection patterns.

Within networks, nodes can belong to different categories, and these categories can potentially affect the propensity for node interactions. For example, nodes can have different sex categories, and the propensity to interact with nodes of the same sex can be higher than with nodes of different sexes. To model the propensity for interaction between nodes based on the categories they belong to, we can use a stochastic block model approach.

Considerations

Caution

We consider predefined groups here, with the goal of evaluating the propensity for interaction between nodes within each group.
In addition to the block model(s) being tested, we need to include a block where all individuals are considered as belonging to the same group (Any in the example). This allows us to assess whether interaction tendencies differ between groups or if the propensity to interact is uniform across all individuals.

Example

Below is an example code snippet demonstrating a Bayesian network model using the stochastic block model approach. The data is identical to the Network model example, with the addition of covariates Any, Merica, and Quantum, representing the block membership of each node. This example is based on Ross, McElreath, and Redhead (2024).

Code

# Setup device------------------------------------------------
from BI import bi, jnp

# Setup device------------------------------------------------
m = bi(platform='cpu', rand_seed = False)
# Simulate data ------------------------------------------------
N = 50
individual_predictor = m.dist.normal(0,1, shape = (N,1), sample = True)

kinship = m.dist.bernoulli(0.3, shape = (N,N), sample = True)
kinship = kinship.at[jnp.diag_indices(N)].set(0)

category = m.dist.categorical(jnp.array([.25,.25,.25,.25]), sample = True, shape  = (N,))
N_grp, N_by_grp = jnp.unique(category, return_counts=True)
N_grp = N_grp.shape[0]

def sim_network(kinship, individual_predictor,category):
  # Intercept
  B_intercept = m.net.block_model(jnp.full((N,),0), 1, N, sample = True)
  B_category = m.net.block_model(category, N_grp, N_by_grp, sample = True)

  # SR
  sr = m.net.sender_receiver(
    individual_predictor, 
    individual_predictor, 
    s_mu = 0.4, r_mu = -0.4, sample = True)

  # D
  DR = m.net.dyadic_effect(kinship, d_sd=2.5, sample = True)



  return m.dist.bernoulli(
    logits = B_intercept + B_category + sr + DR, 
    sample = True
    )


network = sim_network(m.net.mat_to_edgl(kinship), individual_predictor, category)

# Predictive model ------------------------------------------------

m.data_on_model = dict(
    network = network, 
    dyadic_predictors = m.net.mat_to_edgl(kinship),
    focal_individual_predictors = individual_predictor,
    target_individual_predictors = individual_predictor, 
    category = category
)


def model(network, dyadic_predictors, focal_individual_predictors, target_individual_predictors,category):
    N_id = focal_individual_predictors.shape[0]

    # Block ---------------------------------------
    B_intercept = m.net.block_model(jnp.full((N_id,),0), 1, N_id, name = "B_intercept")
    B_category = m.net.block_model(category, N_grp, N_by_grp, name = "B_category")

    ## SR shape =  N individuals---------------------------------------
    sr =  m.net.sender_receiver(
      focal_individual_predictors,
      target_individual_predictors, 
      s_mu = 0.4, r_mu = -0.4
    )

    # Dyadic shape = N dyads--------------------------------------  
    dr = m.net.dyadic_effect(dyadic_predictors, d_sd=2.5) # Diadic effect intercept only 
    m.dist.bernoulli(logits = B_intercept + B_category + sr + dr, obs=network)

m.fit(model, num_samples = 500, num_warmup = 500, num_chains = 1, thinning = 1)
m.summary()

jax.local_device_count 32

  0%|          | 0/1000 [00:00<?, ?it/s]warmup:   0%|          | 1/1000 [00:01<19:50,  1.19s/it, 1 steps of size 2.34e+00. acc. prob=0.00]warmup:   1%|▏         | 14/1000 [00:01<01:07, 14.62it/s, 31 steps of size 8.90e-03. acc. prob=0.64]warmup:   3%|▎         | 26/1000 [00:01<00:34, 28.47it/s, 63 steps of size 1.34e-01. acc. prob=0.74]warmup:   5%|▌         | 51/1000 [00:01<00:15, 63.21it/s, 127 steps of size 5.56e-02. acc. prob=0.76]warmup:   8%|▊         | 76/1000 [00:01<00:09, 97.26it/s, 31 steps of size 1.46e-01. acc. prob=0.77] warmup:  10%|▉         | 98/1000 [00:01<00:07, 122.72it/s, 63 steps of size 9.14e-02. acc. prob=0.77]warmup:  12%|█▏        | 120/1000 [00:01<00:06, 144.83it/s, 31 steps of size 1.30e-01. acc. prob=0.77]warmup:  15%|█▍        | 147/1000 [00:01<00:04, 175.04it/s, 63 steps of size 1.65e-01. acc. prob=0.78]warmup:  17%|█▋        | 171/1000 [00:02<00:04, 191.70it/s, 31 steps of size 1.20e-01. acc. prob=0.77]warmup:  19%|█▉        | 194/1000 [00:02<00:04, 198.03it/s, 63 steps of size 1.48e-01. acc. prob=0.78]warmup:  22%|██▏       | 217/1000 [00:02<00:03, 204.59it/s, 127 steps of size 6.60e-02. acc. prob=0.78]warmup:  24%|██▍       | 240/1000 [00:02<00:03, 210.69it/s, 31 steps of size 1.29e-01. acc. prob=0.78] warmup:  26%|██▋       | 264/1000 [00:02<00:03, 216.52it/s, 127 steps of size 7.61e-02. acc. prob=0.78]warmup:  29%|██▉       | 288/1000 [00:02<00:03, 222.74it/s, 18 steps of size 4.58e-02. acc. prob=0.78] warmup:  31%|███       | 311/1000 [00:02<00:03, 213.64it/s, 127 steps of size 7.67e-02. acc. prob=0.78]warmup:  34%|███▍      | 339/1000 [00:02<00:02, 231.55it/s, 31 steps of size 2.06e-01. acc. prob=0.78] warmup:  37%|███▋      | 369/1000 [00:02<00:02, 250.55it/s, 63 steps of size 8.97e-02. acc. prob=0.78]warmup:  40%|████      | 402/1000 [00:02<00:02, 272.29it/s, 63 steps of size 8.39e-02. acc. prob=0.78]warmup:  43%|████▎     | 433/1000 [00:03<00:02, 282.02it/s, 31 steps of size 1.18e-01. acc. prob=0.79]warmup:  46%|████▌     | 462/1000 [00:03<00:02, 266.14it/s, 127 steps of size 6.96e-02. acc. prob=0.78]warmup:  49%|████▉     | 489/1000 [00:03<00:02, 240.85it/s, 31 steps of size 2.24e-01. acc. prob=0.79] sample:  51%|█████▏    | 514/1000 [00:03<00:02, 232.57it/s, 63 steps of size 8.49e-02. acc. prob=0.93]sample:  54%|█████▍    | 538/1000 [00:03<00:02, 227.02it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  56%|█████▌    | 561/1000 [00:03<00:01, 225.19it/s, 63 steps of size 8.49e-02. acc. prob=0.93]sample:  58%|█████▊    | 584/1000 [00:03<00:01, 222.42it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  61%|██████    | 607/1000 [00:03<00:01, 219.78it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  63%|██████▎   | 630/1000 [00:03<00:01, 216.51it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  65%|██████▌   | 653/1000 [00:04<00:01, 218.68it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  68%|██████▊   | 676/1000 [00:04<00:01, 219.76it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  70%|██████▉   | 699/1000 [00:04<00:01, 219.60it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  72%|███████▏  | 722/1000 [00:04<00:01, 221.24it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  74%|███████▍  | 745/1000 [00:04<00:01, 221.91it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  77%|███████▋  | 768/1000 [00:04<00:01, 220.95it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  79%|███████▉  | 791/1000 [00:04<00:00, 218.50it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  81%|████████▏ | 814/1000 [00:04<00:00, 220.07it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  84%|████████▎ | 837/1000 [00:04<00:00, 218.10it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  86%|████████▌ | 860/1000 [00:05<00:00, 218.42it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  88%|████████▊ | 882/1000 [00:05<00:00, 217.99it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  90%|█████████ | 904/1000 [00:05<00:00, 216.67it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  93%|█████████▎| 926/1000 [00:05<00:00, 216.37it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  95%|█████████▍| 948/1000 [00:05<00:00, 213.90it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  97%|█████████▋| 970/1000 [00:05<00:00, 215.33it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample:  99%|█████████▉| 992/1000 [00:05<00:00, 215.51it/s, 63 steps of size 8.49e-02. acc. prob=0.94]sample: 100%|██████████| 1000/1000 [00:05<00:00, 176.32it/s, 63 steps of size 8.49e-02. acc. prob=0.94]
arviz - WARNING - Shape validation failed: input_shape: (1, 500), minimum_shape: (chains=2, draws=4)
/home/sosa/work/3.12venv/lib/python3.12/site-packages/arviz/stats/diagnostics.py:991: RuntimeWarning:

invalid value encountered in scalar divide

/home/sosa/work/3.12venv/lib/python3.12/site-packages/arviz/stats/diagnostics.py:991: RuntimeWarning:

invalid value encountered in scalar divide

	mean	sd	hdi_5.5%	hdi_94.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
b_B_category[0, 0]	-5.00	2.07	-8.57	-2.06	0.07	0.09	923.76	426.62	NaN
b_B_category[0, 1]	-6.59	2.28	-10.12	-2.80	0.06	0.10	1349.49	453.84	NaN
b_B_category[0, 2]	-6.65	2.38	-10.36	-3.05	0.07	0.12	1179.34	284.27	NaN
b_B_category[0, 3]	-6.22	2.37	-9.77	-2.31	0.07	0.13	1057.71	264.83	NaN
b_B_category[1, 0]	-7.20	2.08	-10.24	-3.82	0.06	0.09	1349.49	395.60	NaN
...	...	...	...	...	...	...	...	...	...
sr_rf[48, 1]	-0.05	1.33	-1.80	2.40	0.04	0.11	1349.49	247.22	NaN
sr_rf[49, 0]	0.03	1.79	-2.79	2.67	0.07	0.10	748.76	384.71	NaN
sr_rf[49, 1]	-1.03	1.20	-2.95	0.75	0.07	0.05	290.28	453.84	NaN
sr_sigma[0]	1.64	0.58	0.75	2.57	0.03	0.02	349.74	428.71	NaN
sr_sigma[1]	1.32	0.65	0.30	2.37	0.05	0.02	178.29	117.33	NaN

5131 rows × 9 columns

![](travaux-routiers.png){fig-align="center"}

# Setup device------------------------------------------------
using BayesianInference

# Setup device------------------------------------------------
m = importBI(platform="cpu", rand_seed = false)

# Simulate data ------------------------------------------------
N = 50
individual_predictor = m.dist.normal(0,1, shape = (N,1), sample = true)

kinship = m.dist.bernoulli(0.3, shape = (N,N), sample = true)
kinship = kinship.at[jnp.diag_indices(N)].set(0)

category = m.dist.categorical(jnp.array([.25,.25,.25,.25]), sample = true, shape  = (N,))
N_grp, N_by_grp = jnp.unique(category, return_counts=true)
N_grp = N_grp.shape[0]


function sim_network(kinship, individual_predictor)
  # Intercept
  alpha = m.dist.normal(0,1, sample = true)

  # SR
  sr = m.net.sender_receiver(individual_predictor, individual_predictor, s_mu = 0.4, r_mu = -0.4, sample = true)

  # D
  DR = m.net.dyadic_effect(kinship, d_sd=2.5, sample = true)

  return m.dist.bernoulli(logits = alpha + sr + DR, sample = true)

end

network = sim_network(m.net.mat_to_edgl(kinship), individual_predictor)


# Predictive model ------------------------------------------------

m.data_on_model = pydict(
    network = network, 
    dyadic_predictors = m.net.mat_to_edgl(kinship),
    focal_individual_predictors = individual_predictor,
    target_individual_predictors = individual_predictor, 
    category = category
)


@BI function model(network, dyadic_predictors, focal_individual_predictors, target_individual_predictors,category) 
    N_id = focal_individual_predictors.shape[0]

    # Block ---------------------------------------
    B_intercept = m.net.block_model(jnp.full((N_id,),0), 1, N_id, name = "B_intercept")
    B_category = m.net.block_model(category, N_grp, N_by_grp, name = "B_category")

    ## SR shape =  N individuals---------------------------------------
    sr =  m.net.sender_receiver(
      focal_individual_predictors,
      target_individual_predictors, 
      s_mu = 0.4, r_mu = -0.4
    )

    # Dyadic shape = N dyads--------------------------------------  
    dr = m.net.dyadic_effect(dyadic_predictors, d_sd=2.5) # Diadic effect intercept only 
    m.dist.bernoulli(logits = B_intercept + B_category + sr + dr, obs=network)
end 

m.fit(model, num_samples = 500, num_warmup = 500, num_chains = 1)
m.summary()

Mathematical Details

Main Formula

The model’s block structure can be represented by the following formula. Note that the sender-receiver and dyadic effects are not represented here, as they are already accounted for in the Network model chapter:

G_{ij} \sim \text{Poisson}(Y_{ij})

\log(Y_{ij}) = B_{k(i), k(j)}

where:

B is a matrix of intercept parameters unique to the interaction of categories. For example, if there are three groups, then B will be a 3x3 matrix where each element give the rate an individual in group k interacting with an individual in group l.
We use the function k, to return the group identity (i.e., the block) of individual i.

Defining formula sub-equations and prior distributions

To account for all link rates between categories, we can define a square matrix B as follows: the off-diagonal elements represent the link rates between categories i and j, while the diagonal elements represent the link rates within category i.

B_{i,j} = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,j} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,j} \\ \vdots & \vdots & \ddots & \vdots \\ a_{i,1} & a_{i,2} & \cdots & a_{i,j} \end{bmatrix}

As we consider the link probability within categories to be higher than the link probabilities between categories, we define different priors for the diagonal and the off-diagonal. Priors should also depend on sample size, N, so that the resultant network density approximates empirical networks. Basic priors could be:

a_{k \rightarrow k} \sim \text{Normal}\left(\text{Logit}\left(\frac{0.1}{\sqrt{N_k}}\right), 1.5\right)

a_{k \rightarrow \tilde{k}} \sim \text{Normal}\left(\text{Logit}\left(\frac{0.01}{0.5 \sqrt{N_k} + 0.5 \sqrt{N_{\tilde{k}}}}\right), 1.5\right)

where:

k \rightarrow k indicates a diagonal element.
k \rightarrow \tilde{k} indicates an off-diagonal element.

Note(s)

Note

By defining this block model within our network model, we are estimating assortativity 🛈 and disassortativity 🛈 for categorical variables.
Similarly, for continuous variables, we can generate a block model that includes all continuous variables.

Reference(s)

Ross, Cody T, Richard McElreath, and Daniel Redhead. 2024. “Modelling Animal Network Data in r Using STRAND.” Journal of Animal Ecology 93 (3): 254–66.